{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Scalable GP Classification (w/ KISS-GP)\n",
    "\n",
    "This example shows how to use a `GridInterpolationVariationalStrategy` module. This classification module is designed for when the function you're modeling has 2-3 dimensional inputs and you don't believe that the output can be additively decomposed.\n",
    "\n",
    "In this example, the function is checkerboard of 1/3x1/3 squares with labels of -1 or 1\n",
    "\n",
    "\n",
    " Here we use KISS-GP (https://arxiv.org/pdf/1503.01057.pdf) to classify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "from math import exp\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We make an nxn grid of training points\n",
    "# In [0,1]x[0,1] spaced every 1/(n-1)\n",
    "n = 30\n",
    "train_x = torch.zeros(int(pow(n, 2)), 2)\n",
    "train_y = torch.zeros(int(pow(n, 2)))\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        train_x[i * n + j][0] = float(i) / (n - 1)\n",
    "        train_x[i * n + j][1] = float(j) / (n - 1)\n",
    "        # True function is checkerboard of 1/3x1/3 squares with labels of -1 or 1\n",
    "        train_y[i * n + j] = pow(-1, int(3 * i / n + int(3 * j / n))) + 1 // 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gpytorch.models import AbstractVariationalGP\n",
    "from gpytorch.variational import CholeskyVariationalDistribution\n",
    "from gpytorch.variational import GridInterpolationVariationalStrategy\n",
    "\n",
    "class GPClassificationModel(AbstractVariationalGP):\n",
    "    def __init__(self, grid_size=10, grid_bounds=[(0, 1), (0, 1)]):\n",
    "        variational_distribution = CholeskyVariationalDistribution(int(pow(grid_size, len(grid_bounds))))\n",
    "        variational_strategy = GridInterpolationVariationalStrategy(self, grid_size, grid_bounds, variational_distribution)\n",
    "        super(GPClassificationModel, self).__init__(variational_strategy)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.RBFKernel(\n",
    "                lengthscale_prior=gpytorch.priors.SmoothedBoxPrior(\n",
    "                    exp(0), exp(3), sigma=0.1, transform=torch.exp\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "        \n",
    "    def forward(self,x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        latent_pred = gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "        return latent_pred\n",
    "\n",
    "\n",
    "model = GPClassificationModel()\n",
    "likelihood = gpytorch.likelihoods.BernoulliLikelihood()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 0.912\n",
      "Iter 2/50 - Loss: 9.193\n",
      "Iter 3/50 - Loss: 4.586\n",
      "Iter 4/50 - Loss: 5.161\n",
      "Iter 5/50 - Loss: 5.233\n",
      "Iter 6/50 - Loss: 3.804\n",
      "Iter 7/50 - Loss: 2.513\n",
      "Iter 8/50 - Loss: 2.269\n",
      "Iter 9/50 - Loss: 2.510\n",
      "Iter 10/50 - Loss: 2.387\n",
      "Iter 11/50 - Loss: 1.988\n",
      "Iter 12/50 - Loss: 1.673\n",
      "Iter 13/50 - Loss: 1.548\n",
      "Iter 14/50 - Loss: 1.498\n",
      "Iter 15/50 - Loss: 1.391\n",
      "Iter 16/50 - Loss: 1.222\n",
      "Iter 17/50 - Loss: 1.086\n",
      "Iter 18/50 - Loss: 1.007\n",
      "Iter 19/50 - Loss: 0.956\n",
      "Iter 20/50 - Loss: 0.911\n",
      "Iter 21/50 - Loss: 0.839\n",
      "Iter 22/50 - Loss: 0.808\n",
      "Iter 23/50 - Loss: 0.788\n",
      "Iter 24/50 - Loss: 0.779\n",
      "Iter 25/50 - Loss: 0.756\n",
      "Iter 26/50 - Loss: 0.739\n",
      "Iter 27/50 - Loss: 0.725\n",
      "Iter 28/50 - Loss: 0.715\n",
      "Iter 29/50 - Loss: 0.690\n",
      "Iter 30/50 - Loss: 0.679\n",
      "Iter 31/50 - Loss: 0.661\n",
      "Iter 32/50 - Loss: 0.644\n",
      "Iter 33/50 - Loss: 0.639\n",
      "Iter 34/50 - Loss: 0.614\n",
      "Iter 35/50 - Loss: 0.597\n",
      "Iter 36/50 - Loss: 0.579\n",
      "Iter 37/50 - Loss: 0.576\n",
      "Iter 38/50 - Loss: 0.559\n",
      "Iter 39/50 - Loss: 0.542\n",
      "Iter 40/50 - Loss: 0.529\n",
      "Iter 41/50 - Loss: 0.521\n",
      "Iter 42/50 - Loss: 0.504\n",
      "Iter 43/50 - Loss: 0.501\n",
      "Iter 44/50 - Loss: 0.483\n",
      "Iter 45/50 - Loss: 0.466\n",
      "Iter 46/50 - Loss: 0.458\n",
      "Iter 47/50 - Loss: 0.454\n",
      "Iter 48/50 - Loss: 0.444\n",
      "Iter 49/50 - Loss: 0.440\n",
      "Iter 50/50 - Loss: 0.427\n",
      "CPU times: user 11.8 s, sys: 21.5 s, total: 33.3 s\n",
      "Wall time: 6.75 s\n"
     ]
    }
   ],
   "source": [
    "from gpytorch.mlls.variational_elbo import VariationalELBO\n",
    "\n",
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "# n_data refers to the amount of training data\n",
    "mll = VariationalELBO(likelihood, model, train_y.numel())\n",
    "\n",
    "def train():\n",
    "    num_iter = 50\n",
    "    for i in range(num_iter):\n",
    "        optimizer.zero_grad()\n",
    "        output = model(train_x)\n",
    "        # Calc loss and backprop gradients\n",
    "        loss = -mll(output, train_y)\n",
    "        loss.backward()\n",
    "        print('Iter %d/%d - Loss: %.3f' % (i + 1, num_iter, loss.item()))\n",
    "        optimizer.step()\n",
    "        \n",
    "# Get clock time\n",
    "%time train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 1.5)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAADDCAYAAACYh0R4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXmQJcd93/nJekf39FzdM5gDg5MzOIiDBDh84IGLXBEgzUuUbIpceh1rxUYQduw6HOFQiKKPWJiOcJjQHxu74dhYU9oNSzJlCyZFgYdEUCJ1WRIpqwGIAAjiIEAQ58xgjp67u997lftHZVZV/jLz1XvT3TPTjPxGdP9eZeXxy8pf/jLzV5m/UlprEhISEsZFdqEZSEhIWF9ISiMhIWEiJKWRkJAwEZLSSEhImAirojR6vd7+EfceMPS+WtjHe73ePb1e7zOrUX5CQsL5w4qVRq/Xuwf40ogo9/V6veeBF0z8/QDz8/PfBhZGKZyEhISLDytWGqbzvzAiyqfn5+f3mXgAnwQWzO8XgHtWykNCQsL5w/mwaewVS5FZ4Gjt/vbzwENCQsIqob3WBczPz/8qQK/Xu9csZRISEtYx1lRpGOPn0fn5+S8DR4C9FEuTbSbKrAkP4rOf/WzarpqQcIHw+c9/XoXC10Rp9Hq92fn5+QVgnsresQ/4ggnrmbC9wLf9HCp87nOfayzv0KFD7Ny585z5PR+42HlM/K0cFzuPk/B3//33R++txtuTjxek9/Fa8HcA5ufnHwU+Ye49Pz8//6gJs29dFux1QkLC+sCKZxpm6fFlEfb22u9fC6TxwhISEtYH0o7QhISEiZCURkJCwkRISiMhIWEiJKWRkJAwEZLSSEhImAhJaSQkJEyEpDQSEhImQlIaCQkJEyEpjYSEhImQlEZCQsJESEojISFhIiSlkZCQMBHOh2Ph+8zfA7Uwz9lwQkLC+sCaOhY2975tTrXurXnucpwNJyQkrB+stWPhvVSOg18w1+A7G05ISFgnWFN3f8Jvxn7gQfPbzjr2Wx+iCQkJ6wNr7lgYSptH6bVLOhseNeM4dOhQY/4LCwuNcS40LnYeE38rx8XO42rxd16UBnDP/Pz8r0DU2XAU4/o0vJh9M1pc7Dwm/laOi53H1eBvTV659nq92drv+2ozi3soHAvbmcU+c52QkLBOsKaOhY2SeKDX6z3f6/WOQdjZ8Ep5SEhIOH9YU8fCxlYxF0iTHAsnJKxTpB2hCQkJEyEpjYSEhImQlEZCQsJESEojISFhIiSlkZCQMBGS0khISJgISWkkJCRMhKQ0EhISJkJSGgkJCRMhKY2EhISJkJRGQkLCRFiVo/G9Xm9/7OCZOZi2QM3hTigsISFhfWCtfYTuh/Lg2kKv19sfClspDwkJCecPa+0j9JMUMwpMnHsiYQkJCesEa+25axY4WrveHgk7Z5w+/X0WFr7Jxo0f4vjxP2Xr1vdx/PgfMjv7QY4d+wPm5j7MsWNfZ27uYxw79lVDH2Lbtp/j6NGvsm3bRzl69OvMzX3ExP8gCwsPs3Xr+zl+/Dts3fpeTpz4b2zefDunTv01Gze+nTNn/pYNG25icfFZpqbexPLyK3Q6uxgMjtJqbSHPz6JUG8jQeoks28TJkz9hamofS0uvMj19NWfPPseGDTdw5syTbNz4Nk6d+hs2b343J0/+BVu2vIfjx/+Y2dl7WVj4lqnL7xtev8a2bR/j6NF6XT7G0aNfY27uoxw79g0umfkwM//xDxh+8IO0vvUt8nvuIfuTPyG/+26yv/gL8ne9CzU/j771VtQTT6BvuIGNjz2Guvlm1GuvoXfsQB0/Dhs3wvIyKAXtNpw5g966FXX4MHr3btTLL6P37UM9/TT65ptR3/8+utdDfe976DvuIPvzPyd/73vJvvMdhu9/P62HH2b44Q/T+sY3GH70o7S+9jWGH/sYrYceYvhzP0frq18twr/+9YJ+4xtF/G9+k3avR+u//Bfy972vqMtdd6G++130O96BeuQR9FvfinrqKfR116F+/GP05Zdz6OaDtNvbyPOTKDUNDNE6J8umGA5P0W7P0e8fotu9jKWlF5mevpazZ59iZuatnD79GJs23cbJk99ly5a7arL1R8zOfoBjx75pZOZrzM39LMeOfQ2t7+DVV3+zbI9Ktop2mZv7MAsL32Tr1g9w/PgfsXXrz3DixJ+xefOdnDz5XTZtuo3Tpx9jZuZmzp59munpfSwvv0SncymDwWFarVny/DRKdQHQuk+WzTAcLtBuX0K//zrd7pUsLj7Phg3Xc+bMD2qy9S4WFh5mx45fNrJ57jhf7v7OGU0+Ql966X+m33+GEye+SL//IzqdG+j3f8jrr/82y8tPcuDAgywvf5+DB7/C0tKjHDz4eywtPcLBgw+xtPQIBw78LsvLj/H661+i3/8+Bw78Z5aXn6TT+U36/adpt69hMPgRrdaVDIcvkWV7yPPXUGo7Wh8BtgLHgRngDNAFloEWoIABSk2j9SJvvLEZrU+SZdvI86O0WjsZDg/Ral3OcPgK7fabGAx+TKdzLf3+c7z22m+Yuvxnlpef4ODBL7G09LccOPAVlpfjdeHbD7L33z7O4ItfpP3UUyxddx3tZ59lee9e2i+8wPIVV9B5+WX6u3fTOXCAwfbtbDtyhOGWLbROnCDfsIHs7Fl0twv9fqE0Wi1Uv08+M0N25kwZd7h9O60jRxjs3k37wAH6Nm9blil7+Fu/Rfuppxj+zu/Qfvxxhl/+Mu3HHmPwe79H+5FHGDz0UEG/8hXajz7K4Hd/t7j/4IO0H3+cLddfT/uZZ8q69N/0Jjo//jH9yy6j8+qrDHbton3wIMO5OVrHjrF0ySZ+9KVTwDSwCHSAIaDN72WU2ojWp8myreT5cbLsEvL8MK3WpQyHr9NuX8lg8BKdzj76/efpdN5Mv/80r7/+n1he/kEpW/a5t1q3MBx+nwMHHmJ5+ZGynQ4c+JKJ9zssLz9Bp/Nb9Ps/rMnWVQyHP6nJ1g60foNifF0ANgKngSlgycgWwLCULaU2ofUpsmyOPD9Gq7WD4fANWq3LGA5fpd2+isHgJ3S7b2Zm5o4V9cm1VhoLwDbze5bCJyiRsCCafBq+9NIZ8+ssAEqdMbS4zjI33F7HqM2noqcB0Pqsoba8ZUP7guaGavMHWg9FnGUTvmR4WwyWbetQ3S/KbrUWR9a1s2jpokPbZ886tGXpUsFHtlzwpfr9gvOh4VtryHO0uQeQWWrSZKYMm1dLlNU5K3gy1115X9IzZ4K0ZamJZ8sv+RmG2sW2iW2jgYhj23TJULdd7PO215XsVLI1HEKr5YbL9vNlzMqWpbZ8KWOWX41FJVsDkWZZlF3QrVtnmJtbmZ/QtfYR+iCV4+C9FL5BQ2GrCOVQrTNDWyNplU7kpjJDlUO1tvG1SCGva3e0pVqE22uZ1jaPcq7z3F5bnW/rYu53lJMKHeZJKVWUaKnWI7gvE7lxDW0qC/PcLNWWtgrelaG02278LAumL6nNvsa7rv2PtWvBqnboqLYr4tn2sDJlZaFtqJSpVjC+LaeqihLUK3kkX6G6VNSNp1SLlWJNfYTa17DmDcvC/Pz8o6GwlfIA9Q4pqRSIXFDZIG6L5Xnu5FcrMcZJ4P54QinL9oVaxFbaoWW6XJQmpbDW8VRRQEF9afUrXrtWggbjyzzN81Q2Xvl8TV3y3I0vG1TCKpUsKxaDG2FxFxy7VSr1UHuEi/KLtnnFZMaVKdkeUkn4+ctBRPIZqktMpqwiI4hY+CRYUx+h5rfnD3Q1fYRK7Vw1WNGAWWYbMtywVbhUKhZ2pmFieQ3qcRQNU0qhdZFXvfFUqLM6PNprtw52alpNUQvaGrodrRRKUQmd5yhzv64AdI3W4yuKGYKliLhBxSNHPhts71slYZWHvbZLI6FEPAVay18Bg47mr79IA+TzdtunDI1OblxlXcmclKFhMH4MzbKVmbyVF8fKlKQ+hqHAifBTuyPUX0aUd0x4k8oNN7TfwaWSic8wmspsEtJ48tFLAm/pEJjau9kF8otM5eKT/whLNl1DmV6+si4i3E8Xeiax0blpJuh3UPBnCv7A4qav4sUGiRg/MVrP2+VF5lnxtvKpxk+h0pCTZtceEL8Oh/uzbdnwsTWi8vKq0oTXrtUa2b2W6X3aduLpTNgPqgJkZQqaCTHIMndpY/50PU3EtuDlLXi1ZdmZirVlqMy1HenasiMIqcRKpRJrf789IgusEYjJxOgBqnpEsQHLVQbWjubzW8/HldeYfUTK0mrYNC76V67NCI/CVcNII5Rr+PQ7qruO9DW1S8cdnUJl+DaL8KylipcJag1vrgBp0c9kh61zF6QNo5ET17xVqRtGCVEBJZSBFsokaoeJ8WP5GDmLiN2TrMrn7nHv0MrwaeskDNO5G99X+lI5NM18/Do02/KsQroIDKEXGrFZbtUg0oYxFNfu1LKaDVi7QXGd57KDTzLFjTW+2xFi/S1Wl8o+Y+tk4pdmmkiGZv1fH68cWn+oWld/tbR146kyPI4cr61tomZPKS61c1/FeJYG0hp/4fLPZXniFukXJWVFylAuqBtf5mNlSsqWnGGGVfv4g1W9TpX969yx7pWGXB9Ko5QcxbWWVGpiVzj95UhsVJD3nVxMHuFwfwQyOWWx6axVEkNx34QPRQf1ihX5Bqb6IduBEtdOXjJcQDeVKd+iRJZSWuYvZywTLD1iy0W5nJDwZ4y2HVzZqpS6zdetkz9QxZTEOHWRPIfrlpQGUHV2OWq7Dzy2Dp3c2CgjjBpfw2kqXptGPHekik+T3euy40h7gpe8YQkgZxwjMJlVYIz4kSmkqim2MF9N1/69WIfzzTfuLLRMFe2wlkU5g5DKYLyZ6DgYsYoz91fe5X8KlIZ88K5BzY8nN+jI9ah8JFawpMEpxk9dIMKGrLgFO8azpHITkeC93cCkLD9m68gySkNozTjqpG6YzZSKS2zS0kKhlQZRe99u9mqqgzfziDWMvxgb3z4l82ySISlrrh2h6rgx2Yp1y7psjV7axge3yRWRxLpXGtJYJWcYFcIC06z1rWC4m7ya31z5685zfd1VJRvd4CVvyg3wlgSBhME7wpYxVhrJjLyOvCUpFVLM8mzv2x+xmUbjTNGPFJOhMW26tXxK5iLxw7IX3yEc5reBi4bwNNMILDMK+LYKabNwjVX+Wi82lUSEj9OgUigs7+OmC/MuNw+VG9kGrk0jujyppVYmvqNa672lJsj1ONG3HbGaxDZxiZ6phqI9YkstOdMZOWtoWGJ5bzVseKwTu+1in39lc3IN1f6O0qYZxSgFEL4njef+DCPZNAK2CnedGjOM+m8gZPzR5cTRbKyK2yjCecU3/bjhZd1iyyF5jsPGF282RnMjwkYYU0PhXh6xtySRfL1R2Hbw6HSgeTru25jCz8HfCxGWrRElievwYOIjNnDFeYu9cvVf60+Oda80YvCXKXKkClN/g06Zo6AWsUeoAr9jI1kkCyEs1QYdWWZkvSpH51jHim2kqm3sYtzNXTYPUYa3WIzNHOTmLvl2zP6QyqR85SSfRf0Zjm6H5nW/VA6Ct0id4gZQRHhYZn0Z8+siB5j44bdk08DvxG7Him+IcjfgxA2gUsBkQ46//oydQIzv/QgLRrVZSJ6inPCUK0JctTwpGqlHPa6e8JSr2MTVeMq1idrsS34cLhnVSZpPucqZgNv2459ylQOSOxsu6xDt6M3T3LhsuXlcFJu7er3ex3u93j29Xu8zgXv7e72e7vV6z5u/L5jwBwy9b6XlW8Rmyc2nXIlc2w466SnX0P2x1zZuDl4d3Pxig/3Yp1zNyVCg3KgVZSJwHRyXI8sHufwY+5SrzHfELEkBnTMwdQBmHxu15GhagtqOF+7sfvqmk9MFRr3RHsWHP2vSxHmxecilk1y2nDtWpDTGcBK8bX5+Xs3Pz+8DfgF4wITf1+v1nifuW3RsrN4pVwvZGJkTf/LXc1VYbCTxl0w2XAvqCqXcVHQup1zrtM5BSBFYxVJXMPXTrsHKCUOlJ+qxU66SynQyf8vjkubdn4Kb/o2MWIe3WBKsh5eP0rYUO+VatdNKT7lKPmNnUUIz4nBdLgabxkgnwUaZWPTm5+etkvj0/Pz8PnF/VRHriNVUcTyN6z98md/4p1ybZhxN5oKJDWVWgBoKGmuVG5nKeTs3G1C2yogZTJCnprrIdEG2xmuHeDq3A8qZYNzA7Q44zadcY/zG+Y/5XvGN7xd4psGYToKNw53/WgvaG1vSrBxhO8C5n3KVjWAFRNoVQnyMPlXpKwVbpnvtz1QkdXmInnKNTf0DBlBn5jCJIbSerkaVuI6ecpUH10LGWfCVmOSnTFZ/hrFTx5NCzgjlocjR7RbvuFa2YnJjkQXKsLzIsnB4W0+nXO+tzypqH026t9fr3TNqxtHkWHg4HBhaTAUHg+K637d06NCB2cMgw238waDwx1jZMnLnejiUZ1bkOja+hq4OKsk8h8G62OuqTrlDl5e1oaIueZVPCxjmOS1gkOd0TbiiWApow0eGO1pVI7YuJNIztPlLHA3kwyEZMBwMCmquB4MBXWAwHBbU8NQfDJgC+sOhQweW98GgrAP1Ohiec+tIaDgsVvqWn2B7yHV+rF3s83fpYCBlzG2PwUA78ap2cmVMtn9FLT9NMhWqC04etg5VXYqyFxZO0O+P7lNNWKnSiDkOlihtHcb4edR4/DpC5S80iGbHwi2GQ8iyFnkOrVZBO502S0vQ6RT3K5oxGEC32+LsWff+0hK0222Wl4vDYsNhpfUrr1tKTD8nH62qvDK0hizLTB0ypy5ZVvDWbrdL3mVdhkPodhVnz0K322FxETrGv0bLjNJtQVvGaJgJWtYmtHSwMwHEGGdsHPW86mWXPJgZheSlDDe0I8Jb5m1Ky4ycbUPtdSbKt/xU5oNmA6jVibY9/HaxspU57dLtZsXzLttD0e8X14PBKNnKhGxlQrZiM4zR9QjVRcrW7OyWC+5YOOgkuOZYmF6vJ5XCPJUz4X3m+pwRn2JKrd10ytXmJ41ashxJpT1h1CnXGK9hI9akp1xbZ40hdDFyylUuEWJvJFT4mLsScZy8IuEW5VONveaKvUURyw9vB2p0+TKqPWzU2DICJzyWrqrC6FOuMUOov9SVMiVlKz5AxfeaSNm7wIbQEU6CvyOiviDSfMI4In5+tRwLywav+kHmXMetzFbJxIQtpt1HzTTChqtx19KxzZUxBbXrm4orfxO2fy9iT5CIvYodxUwE45bRLP6jy2885VreLwNGZOkODJJL33wT7tTxAUWWJ9PH9vog7o+PZtE6VztOhdVwLBxyHFx3LPwC8I+a0qwcbkPEGtLfoBM7oejmWz+JKJb5ET7sRH7o5Tm5BVvWRRreiuvOYsaVvwHD68O7KWsMuNexU67WCGri6DzHMZLW8nLsIIGyVW0/BdQMoXbZs9JTrjWqgWwIDKF1CoZbIbCwopIJWXsJyYNsB3eDoO+uz+38o0651pcWzfxU92Obufw0K1ca635HaNMJxdgUr3kDjcVqnnJtShPJSaYzM0x7Dmr2bwu64RWRQL6GG7GbKKz/tF/hGh1ZnaZXqpF9Hk0H4LxlUISfbBne9Qm47dMyZciIOFqGmjtkOJ4/MwzL3vinXMcRICn/Mo+Vd/l17yO01Vf0gdZAMwBa/eJ7VC3zZqG1VND2kmYJaJv1vl33txc1S8rYAzLIFoEMWsvQV0X6QQuyvmKYQSuHgQKl82LZrHUhCyUNcVncU1qjFai8WHJnA8izgg5VwftA1cq2dVi2dSi+/bXnD3Jeux6u/m+aH98FVz6ouf7PQF8nTowayNedssNprcu3J0rQMl4tT1tNp7pN82L5tsXm13TKVcwovLrU7DayDtNHIW/D3F9BexnyrHj+R+90K6G0aY/ctMNAkWfQKttFF+3S10a2iqJbS+7u1izLTXWkS8nYLmS3AzfvPMUPF3WwstUaKgYZZH0YGooCla/8lOu6Vxr7/gMcvgp2PKt44wbY+UPNwZtg9w/gwM1w6ZOa194Ce54o6GWP57x6C1z2t5pXb4HLH8t55Va4/DF49W2wZ17z2tthzzy8vh92PQIHb4Udj8Phm2H7k3D0Bph9VnF8L2z6CZy5DKYPKZZnoX1aMewWjac05B1oL8LyZphegLM7YOYgnLoMtrwMx6+G2Rc1x66B7c/B4ethxzNw6EbY9ZTi4E1w6Q/gtZthzxM5r70Vdv++5sp/C8MbNLt+vaBKB0RMvlLsquL15yZN943iw0Ld4zDcoMlOw3BK014C3VWos6BbRT1Qxe8sL+qj+pB3obUIwylon4LhdKF4hjPalFF81XawWdMF+ltV8Up1TtEClmc1G4DlbRR0zlzP5WwAFrdpZoDlbSZ8VtN9CfpbNFPA8ibNNDDYrOgeLXgPm6AVt/xLCmfLps++/CkYbIDuKVjcBjOH4fRu2PQanLgSZn+iWdgLc8/Dketgx7MY2VIcvLFoj9ffApc+WbQL1+QsXgXdw5r+VugegcWN0D2iXXpMsbwRphY0/Y3QOaEZzJjntwHaZzX9afNcp4uOnneAvjafolXFtiCrezLY92uwvAmmjsPiDpg5pDi1Bza/ojlxFcy+CMf2wfYfFbK1ZSmH966kx/0UKI3t34UdD0J+GVzyEORXwfaHYLhXs+0hyK/PmHsI8hsUcw/B8MaMm74KwxsVN30Nhm9W3PQNGF6nuPH3YbhPcePDMLwKZr8F+R7Y9oeQ74RLHoZ8DnY8DHoT7Pp90DMZ6ozpaMug26r4rKbdYJUDHYXqg55WqEXQmxTqFOitsPs45Jcodv0e5Lthx0OQX27qcrVm+0OQX0NRh+sVc1+F/MbwGwPpP3N5W8aGH8OZqxSbX4Y33qk5+fNw8jbY9Cicvlkx80M4c03GhhdhcU9G9yD0t0H7BOQbFKpv8m4rskUYblK0F2CwQ9E9AEuXK6ZfhMVrYcPTcPZmxczjcHo/bPobOHU7bPoLOPUezeY/hZPvU2z+Npz8AGx+GE5+CDa/BU5+BDZ/A058RLHlG3DiQxlbvgnHP6jY+jAc/RnF3J/Aif8BNv8ZnLxDs/l7cHK/ZtP3ob9LcfMvQ3+LonsGZyZS2meyDPKcK76kUAPQUwq1BHpGFW24RXHpCci3mfbYCTsfKmTgkocgv9K0x16KdriukKnF2xVbLoOtryqOXwaXHFQc3gWXvAaH98D2VxRHLoPtL2uOXg5zP1EcuwK2vQTHLoetL2tOXAabXoHTl8LMAcXZHTB9OGNpFronMvoz0DpbKGs1NG2SaS7570K2NoI6DXqr4tLjkG+nqMuuQraWb1djLXJGYd0rDQ/CwFaumaXBzRraIsfCy+4Xm3Zb93faPe0Z23UZjDtmXUrjoTQi2hOhsg7m/sJ+zTP/Eyy9RdN9EpZuhMUutFqKk3uK6fOJq4vanrjWlHmTw8AI1jT6LZBlmoUbiqIXrikOqS5cDZ0OHL+82Ltw/O/B1JTi+McN/QWYmso4/gmYnm45dGqqzSufKPaenPiFIv4rHy/yf/XvFfmd/HlotzNO/qypy+UFP4/8P0Ab9n8aV1mUjVAzHIbaI9Yusde85rqzkHHFX8HSDS2u+CvI39Liij+H4U0trvjLYmC64rswvEZx+V9DfjVsnIf8csVljxTKadOjkG+DLY+B3gJbHgVmdPHN6I42a26KlY8V1RzoTChb7ZV3+fWvNEYY9wDPOUvokBYQPxgl9gtUb/J03IRh40tjndn/4NkDxjBUBuvSYFnNM1h4W9GBz94KLSt0JZpe+VmjoRJhcZabw8M7aqt+KT2qyXIlrd5knHwzdBfMs41tNw8xh2unGWHRHH0t7TURRPeqBPgpLs6tLqN4XAnW/dsTD1Ko5GlJ+eCbPjjsuci3t13DYElD+xIim6U8oZAbomJ1kUbEiOs8nUnRsdZ7t8M2b/ip5+MaJJv3nLhGPX9iZzua+z1aG9//pITIXbv513mNvmRUKthmZQ5NhmOpHGIyZdoj7lUsXKmobNVnOLU0jozFzrfY+NKN4jlg/SuNiNBGhXrMd6YhQavT8n5M6YTKkm81RnJQq0NTfDsT8RRWLIV89VdmNIKbWOe01D5vRlL5WlNuris5bPyGrhve9DpUFB7Oqymd3Bho0407EykLipQkZ7OBga6cCWmxL0a+fYpUZTWw/pVGrEFkJxenKKO7JeUpy0jHL++3qlODDifG8KbNb0vrZXui27Cr0nP7L20ekXyqR+QqOLtNvTr5GFJPqvYXiuMqCzedzKse33Z2u6nL3SjlU5mvhZ3BVMsUt7jas6+1vdMekyIiM9LmIfeiRHexCnjfsa1tqNM2PPZJiYhRnIC8nivWv00jBjE6eEetzcPzGtR6lrL5iJ2GUK19tYkfFYGAcGgTLmm97Fg+UuERM4Ra1u0msFLJVDsPi2xjthEtfo/qXDYPd8oesz1UutfdLWmv/WPhbh2qa5fn8ryHFudb68o+1h6SynaJQRjZZTvoLDwmh/bNOIOL/D5uyezomZJTh9g3dldBaaz/mYaEVBYRQ2j026EjZhgqRCP2CrlMcdJgxm6ZNraUsnUpj32b23bNLNaptj/qckgId+zYqB2ojPfb33k7Oq+4F3jp7l/Gs+nHVEqWEzl9rzMba1P8edWo9FK2VEy2JuSjlAtLQ8udiAG0zCN2cnkVbBrrf6YR224cM06NaekuhdEGyFE+MNVUiGmhacCoYa1pehwzakkXeUJoD9+p6WyGw3eWjNhK1Dit59swogZZc10FNB3aigy6nlFWnhSVJ0H9k6mZ039Co3g5K6zzT7xdtDAqejzLgIiLQh1RHt5Tjxt+3DqIutjfKpSXvLYNcDG8PRnlWNjc95wIN6U5J4TeWiBFNx7fM1JKIY097IaOP3Ka27S2tUXEjFuRspe3wxO/CqevtCHygJTbEWsZjsHN6PuV/LstEDdUyhlLJHdvZhN5XdxkhAzdayrcQsqYVAaxt2CivPLJNNk4Au3rz/nGhJxJrwArUhpjOBYG4UR4zDSTQz6UiBGqVAaRdain9a2yiA2V+OvhMn3MkDWmNb0MjbjAk3Up19Cl8VDaAWLFxuqmxJ/DlcHoJU1Vhv23yYoZAAAfl0lEQVQhywq7+6+Ms2EXBdVyxsbHua5zZu0F9T/HhjCmgbKsSYMtY2zZsuVJ2YqdOq7LVqyeZWBExi600qDBsbCBdCI8Tprx0aBBpVb3vGeXESOjQ+TjP8E8JF8jRr3JFwR+t9WiTmU8O6J5yw9pJ2jiRgfuhW0aPpUzjFgtRO7RzwaEZ0Wh5WWwNiGDaI2ObI8m2YosP7z0MdkLyMm5yMfItDEFdQ5YqU1jHMfCe42Tnv3GN+hYzogtmnyE7ja+Iq0vSs8nZb9f0Np1K0CHy8vFtfVJWfNRWc/X+tPMzf1cFwelLJVTUftbBeIOtS7yqPm7zKj8X5b+NAN16NZ5FnTz9wZsuAk2HR3wxttgaD5pYG1gcT+no3xThhWQrzxcP5vSR2Xlg9X6MV02tC/osogf9tVa+dUsaK4Ln6HlM9XC36mho9pD18oaKVuCDo3v0LI9rGxZWbOyVZMxR7YsP8Z/ar0O9gSvU49RdYGgbB07coTlhj7VhDU3hEonwpOmb/IR2rK+JCM+KUsfkzbcaPuYj8oynaGZvDbpLPXWp3WqVPXlr1ocO33NzHVJY341RR1KKv1rGrr59Yx3/iK88r+2eONtxfmMwrelKv1FGn86YpBrtASV93xDKI5B1JbRahW+KVst6++0Za4z+v3iut8v7tvrwcClBc9K+Om0fjaFX83AMlTZ148Vk24c0Q5KtHHsebdle1nZMu3l+UUVMunJlnnCmeSvNktQxoGyrbG3d8fybvMSdZnbsgXd0KeasNK5ykjHwr1e7z7j1g8qJ8LjOiNeEapXfOZh2huRb4R6U3xBPUirdMRg4KyZY+vJWHiTcbbBfoPnUj+cfWxvhPsUZF7hjud/xzRYdO1+7AmHy6ulNOW5BWSZfAaB6fq5vs2Sb9REnb0NVTZZ5K2INwOK8SNtHMIuMxbvgueVYKVKo8mxcMiJcDDNOSNi8PFOucrva8hTrlK5xGwltUayxikNzmaaIJv1uFKQm6z4sRO7YqOarEtWHkXR4rZrPIx3XIirzfAyxTeAumXGXORZQ6jdIVq9rAorsrqX+IKa8nPxrEe90ZDtIW1BkfYoc4oYREvZki4LpYzFoP0NaqN070jZquUJrMop1xUpjSbHwiEnwiPSrA6kNo/sy/BOuYpRJDo6iHVxCdlLxLrTxhk59scUVGQjmhZ1k0KuRX/x3xpHXlu6uQSvY/3KD1ciXNpT3BlDRWU+0oZiqftEtf0Vab8g04zRket5SsTaZdQAEkpfKycoW6OMufXooULjQjAxzodj4dD91Xcs7E2L5QhY3nDpqANnEDzlqky+dXEtwwPLHLkpR+YVQ9j0GOA9ch0/5WqvV3bKtW7TkPdjywl/x7ttJ/sKtcntv8jV6zTVtdMuIpNQm9WvHcUTebPWeMo1tutYIiC7QdkyToSkEnFkTNTFQzrlSmhIMsGRqWCsASOvLWPrTW8ZE7qOCZUtI1K2vF/SBgPB+Kdcg+ww3ozDU2FlaQUPlpfRMwd/W7keGX7Op1zr6cZddjQ95wj1yoldN9gVPD7k7NbMNJxcYvYSL/MLb9O48GjS3qKTR0+51izUEF+eVKO4bwCVo9nEp1xjkLYMYbfx8osaI90SyzdAKzrlGiwyGq+K79o4fNtGJqjM18LOYGR+IoZpD8us0x4+cyKDhk4fMWB636WN5SeXlYHNYI4KXckp11VQGuv/7EkMUgnIebE85WohT7laiIYrX6U2GZ5kcD2NMMA1jnBN1nWxEc1+4iD2gagLe8rVXkulMZ5SWpVTroE2DLZLLX6IKa+Tx5RF7O2XVAArOeUakC2bJ5BOuY4DearVO+UaTegLjKJQFqp27RmtQg0aMJ4qCuFx0jaMSN4pV5tMrlPL0QYnpv/ZwHEhlUhdGZSFjkhTX3a496vt4GEq0zcpJRRu+7iZuBRfRcl5VRARg6dc0qqY0TFgg1LmOjivCymvEQbQ2NzQKXMFWP9KI/ZuvPwRsSvEfIWWa2Tt5iNHj4AtRdXi2XuOYonNFGIQvJQQ7uRi3z/NBob1vk3vjur+smXsBVNg9jJ6+TL+KdewbSN+ylU+Q/8Zh57yqPYY+4NNFtLtojx9LJe2Np/YjDEwQwkNUMEn3TQbvRjenlw0CLy1qNOm+I2G1NjSYZy9/LG0DcuRciYhaImGV4DbvwtXDaA7lfHcp6BzWjOcrnJWSrLQvAypFx1iv5JVl+t4VcVMIcKLPztyn3s2hLwF7TMNS4px74UQM5Y32ZRkNtpVhNFlUKB9/TnfmJi0riOw/mcaFpG3HNoayARtOokop/aeINhLu2ZUYg0dNIS6cSpH2uHOWobKEVBs7vJOudrtyUuKN/0m7P4W3Pyv4Nr/s7i95fmCbnrJ0NeKdNMHgRw6x4ABtE8psiVFtqhon85gAN3jCobQPYxJU9Rlw6sFnTF5KnPepVIdtoMIZUzYCDvzasHT1icMfbKgm39Y5Lv52UKJbH6hiD/7FNzyS3DTv65yduZP1hCqVK091OhOWMqSCDbtpVuu7MgNhJ5s2XQ2vrC7lY+oHIiEraUui9IA79Ul3B+Suz/gaE+z7RV4407Fzv8Kh+6GXS/CwZ+BS38EB++BPU/DwQ9k7HkKDn5AsedJeP3vwGVPwoEPwGVPwesfVFz+THG950dw4F7Y8//Bwfdl7PptOHS3ZsdDcPjdcMl34NjbYfa7cOKtsPkHcOZqmH6l8GXROQF5V6Nb0DoD/W2aqQNw9mrNzPNw6ibY/HiRduv3iry2fdvk/dXiq2UZYL9NbT1w6bYVVnNd9sPIrKlcimku+UvId2ve+fehnVGc6wCGGlq5Lg44afNZDbOrUqGN8CmU0mgNKtPkQEuZNJk5HNWCYQ6tTpH3S/8UXr8NLvkevPYO2PGn8Mq7YecfK16+E3b9UcZP7obdfwgv3g27vqX4yXtg98OKF98LV/6O4tr/Dt0ZzdWnodWGYb/4HMNgGTqZoj+AttIMc2hraB0CNkZGYGsH0Bq5HDt+vWbT03DiLZqt87DwNs3sn8PR2zTbvwmHb1fs+DK8cTfs/E9w8L2K3T+Cg+8zsvV+jGzBVU/XZOuDNdl6Gg58QLHnOThwr2bP/1uk3/1FOHi3ZtdX4I13ay75w0Km5/4Sjt+i2PI4nLoGNr4Ii7sovia3AZSGbBEGc9A9BGev0Mz8GE7foNn0pJGtv4Zj+3O2/TEceadm+9fh7GWK6Yl7mYt1rzSe/SXF8v8GnSnND++DzhQ8/Q+hO6155pMwtUHz7M/D9EzOsx+BDRtznv1QQZ/7IEzPaJ77CExtyPnRR4t0z/1d6EwpnvsUtDqaZ36x6AxP/2NQLdC/RNGzciDTLnX6r7a9DzSoDHQOWQvyoekIgyLv4S9DuwtP/RO46tuw87fh4KcUu34DDvx9zaX/AV7/BzmX/nt47R9o9vxfMLgqZ8sToK3xMLIztG6t3/A65JdA5zDkWzWd45Bv1HROg57SdJZAtzRqSLE5LDf1sL/NZxn1NHQWQW+CzinIZ6GzAPkO6LwB1/57uO4A6Cvg2n8B+k2w719Bfk3O3v8d8utzrv7XkN+Qc5WhV38O9Jtzrvo3oK/PyV6G4bUw8wosXwkzr0J+ZUbnNch3w4YDkG/TdI5CvkV8h9Y+9nLmB4fvKLRiPgPtBVjeDdMvw+P/h2iXjmLYL9pj8M8KmXrqH0N3Gn74i4WsPPMpI1t/t5ChZz8KnamcH/9sIWulbH0IpmeG/OgjRrZ+vsjvuf8R2l3Fs/+waP9n/hGolkL/s4BM1WVLU2gMar9rsmVpJVuK4a+YuvxTePPuYVIaWuPY95TShuJQyCEza+MMio5QhdvrKr57LcsBJa6zKryEqhrYmrMz0CatLtNqJ6/XPgw/eT90u5qXfgampuDlO2F6WvPyO2F6OueV22DusObSU6Cu0Ox4Agbd4rup3p6Sple5cmYyYrNZ48q4bhjWphNrc2Rbm8egi3ghauPFtsh7337xDKQuH7qT8dIncoYbim+4Busl2gUQVDvUyoKVtUq2ZLzcyaeSSRVs9/i1rmTLs4QqR7ZQUrYkTYbQAOQbB2ltj+0HkAa7UUbB+v2YWShu57ZGRN9oKOGZ0Ez8osyzlyueuh82DRXHN8GWRcXOx/HWsY0fNIp4M7ffPpV2AMahQc5r17E9C7E3GE3KwZZh07eLb+sOpjUv/i9lhahPB+WxersEa+beFh2WIVm3ak9KTImHZbWy90gZKw1qgbKlbI1Xl0mwYkPoGD5C7zN/D9TCPL+h547w2Fd9R8M2WGynoWuIk68hfXdyOOFV+dKBzQiOxR6HuINfVzFZnvPcDR9szHjlEzDY5BpzS/GIbEpSAYOgQ+vSV+ff3KvTsc9XxE4fy5OhMp2E3NkrNmgNO5q/+XV45P+udyqXt2bHyk3t4spW1U4tES6VgcN6gAZbYyT8fSs4tCrzAp9ybfL3aU6yftscUNtbc8Lj+A1dDfjbMWIPPqblC1QyGh4d/Pzd+OF83bRxyJHJljXKo1a9FHndDGUKio4/sdeBxLqCLEDclZvqxIa1KGL7bQL7Mexs/fQ1sHRZ6Fk1tZ0tQrZDrKbh9onpOzkDsd7HZLn+DKSe//gKpcjT0vHij8JKZxpN/j731sJeoPKjIf2GnjPkyO/vOAx3NH/nYaQjlh/xkSXHHn5IUsKzFVmH6loL6ob7s5uCZkOx3pcb1MS+gJCbgKAisPsKTB7OMiC2R0G+opYbm6Tbf7tByrjMiw6Vgidvw1T5OnuUKhtvih6XLfl5BSlLlrrfpa0GD7f94ofupLKoL60jy8DY8rBs+kHw/iRYUx+h4gj8fgoHPOD7DV0BQlrYf3hSmTRrXKl03DXjeDOXcaeZTUsgxuK5v6Wggw1mNSyXKWVx4tnYH6N2C0ZG+mpiFrY1eG9yZDoZX8SLGmelrcNLV0asFyJoGI2G47KTW1mSqlamd5+rnM3Ww92iY7IV50/y7uvwlds0zosh1CxbSgc80m/oqBlHk2Ph4dA6rnUdzg4GRXi/P3Suq/uS2nha5DtwruVUMg47GlTnQmJLKOks1+cpVkd7Xdx/5d4hZwew4dSQKx+vnCMP8rx0fJtRc2SrVElbNerVsGYI1VTWAUXxtiPojNdSrcsyW1TOkvvDIVOGt1admvz6ec5ULV1u6lg6W7Z1qjnldergzcpse/jtV7XD6Hbxn7v7/Kt20+I6HK9yuuxexxGbgdTr4lLJ+/HjxxkMLqxj4XH9fd4zPz//K1AaP4/Oz89/mcpvaBRNjoVfeilznOVaB7aVI9uOue4yGBS0cGBbXGeZcdrbapcOb4t8OuS56/i27uC2eWVXb9C4IaywNRbhdWe8bl0sbx1Th7ZxuuvWSW3o8vqHYe83TZ2Eg9vSyXKrhTY10EBmFIKlIyf2No5Ng6tIrFPc0hmvuW53OgU1PHTMdafbLcIttfEMbRmaCSfL7Vo5TvmWH2+aHm8P+xalelmkRLtUTo3rsmXbodWaMrRj2qXjXdfbsd1usbxcXFsZK2Sx5ZTv89uM+tuTok7Kkdu5ue3Mzl5Yx8JNPkLp9Xr31WYW9xD2G7pqkK9a/XWnex07I2SnlP7o02REC81EYtNKufyQNGzT8NfWbt1GZFhSO+5KdRYUz7oEa3HK1+ajGlz+2Acpp/aTerkKGGad8kfyEVsu2mXGqAr4j7FqF/k9Wilr4XyqzzIIG1PjZNbO+UKQtZfG9Qt8yrXJR6gJf6DX6z3f6/WO1dI4fkNXwkNo9C6obIg8SKXlO8tcQZJrQf9Vq2y8+iOVCkyOcKImYv1ZOZZxDWdSWci3K/bchzR4enshSjb9jho354aZHfdkaON3aWPKQxp1LWIHBsu+EVchsROz/oldC1dGxpWtypWhm09zeR7HI+oyOm5VxkWwuWuUj1Bjq5gbJ825Im40jKVwG8wfbeJLiTCqw0VxuEJmp47jjmw1LkR4rIO6SkLJY9FNeynGY2Y0Ikok2j3Em51Y+arSqOH4pYIsA6IsyvaQcWOG6aa9PDGbg78NwI3n8xG3XTTXxYsheDh3rHhz18UD2eBSPOWeWisA7qavWEeMfzM01sDKK0ta3Wu5B8uU961CqzauuXUqr9vCdX5ViXD2sW+HGgMoSvkndssEkeWEZ8UPz3LKk7pCwckTvZ6SsK9oI698/Zr67eEr4XERlikrS1X7yA2FMp2UHZz0o2VL1tCWEVNMNu8LvLnrYoDvTs6lFr5/h9h9CbEGb1hy1zjDn8WMK50xoW4KN6Upl0lvGeKxOsJl4Qi7SFCJNEHOfuR+j5iSsCzZH9K1oXJ3tfo8+e0h20XakuTzbWqHioZlyx9oQjzW5aRpMPHTllc6HD5ZnmGse6XhfSynhNy0JY1U0ngo48dmBbH8Q/HOdSoYFm65eSjLwtdqIEb58ocQmJrNQ5lrh9bjBewdwbgN0MJ24RlCBfVmLNKGUTOoOvwElyfjtYe/3MC5ljJjr6v2GDrh1lDqj/5ySR1eMvt8hxSgzcPlVfK+GjaNda80ZJ/2d1GWMU38sDKpryvd/MLlxKaHo+I0GdhikHWQrvHkCOW/RGnIX96PGBfHUg6RZUrJm4wnDZ0jZjZA1D2jZ9Mq33LV6xJbeloaq6GckVhZaWqP8KYuf2cvDdejZh6xZbPIqcw6ffckukyPTfHGP+XaWHIw/3C45EHS0cIq82mqg2TJy19ey/fOobctquGUayxvwaHHU6zukXI8F3lyCVYuf0YtBaSiaVo2hGce459ylSzHBg2Z36hTrrGyCVJZl5Vg3SsNOYpW8t50yjVsAPUt40SofPhrecrVtcPIU67SEKo7ykkd04T1ZUtwQqx1sFPbe3XaVFYJaei0p1ztN0btpyVqvLmFayefUklo4Rovj3XEelajlpb18KZ2aTrlKg2bTbLVtDzxeRz/lGv6hEEJ32gVe/Cjp4bytVnspcD5PeXaJNSinGH4rl9cNWpHx58Rb0Xi3TFQhpghVEtsMbWPLcolH2KzGFnmbjarEgoqf8fRZFyPLXFjS105Exn/lKssTwfujcb4BvxmrHulIUd+fwONe12Fy5OJMp650nKUCMercRQNa57R23haUMS13JgmT7uKOkh25Kh0Pk+5yjKbaJOUeyOsdviqEF8uxtRffKofaxepRKSx3eWxyic8y/X5DJ1ylTyPrstqnHJd90rDX4+6DRFTHuO9Mq03eHhdO1qZTKZofCUR7ngxXmOoL0McWgab67U45RrhxX+MbvxJT7nK/GtarBYw3nKk+ZSrSOVr5ZH5xb7VMr5sxeviKyTJc7Jp0LTe9OloI6JvC3HLqSzncrOYhArck2VFktoShTFXGkDttQzXLaFsqgxdGvlUZd3GYY2gUUPomDMBqVzsTCDm7t/br9GQv8eXl6zeHu4ywYvphY+WMbnciMlafHOffMsSk9F6+WFZihlh/bzPHeteacTWatUBM/vwZCfPxH3MtXtOoDqLEmv42Ojsrzv9zUOSyrooES7X2HJnaHg0Dvnh1EVlg0ZN52EGlIIc42ReQcRmN5HvzoyhUV1+Yt8/HWEHiLEaaxeZrqLS4DnayB7/Ol2I51C4H+b7Ic3FfTnYnTvWvdKoIKeCxUO0nT7mZck/5eouayqlYq+b1jXxho2PWOauJ4xyqhleO1d1ksWGlxTW8Fk/mapq4T77rhKRcWVeQXhvuUx7TXrKNZCvU/5IPoLqjkp2ZPzR7SM3b8W/S+uWLwcoeZJ6vKVzTKGEZ9Crecp1xRvRzWnVBSJeuEL3m9JMhtg0TY4W9runTadcca6rDlqdkSiyjLVshu/mrRpRirTh5UnstW4VT25Ei5xyzV3pdpYaNaotJwEpVYEaOmEyL2EglSi7Z+wtyApPuZZ1KfOtcx2G3x6uH41YLarZp2Ut3A7+KVd3QFrtU67uo43M7C70jtAxHAt795vSTArfaBhey/np3Ab0Rxv50EcpiXo+IcgpogmNOqqNQdYt0kHFOt8zdI65hBCZTsAnnlIpebG3JynbiSaWNbG1gxjdQ/C/ZB9mKb5ZS3b+MEtVu42eafpv+qSMxhGzF/vxJmzHANbasXDoflOaCREbtbMItZMrdyOObyCNGR3d66pBQycTw3lUPFoaMnTV4RpnfSPshKdcZbg85Vo3eto/77u0kyk82dmlIVSJ76B6hlBZnlWA8r7Nryq5RjNDw0oAL0043JeB2ClXadOQMuYqLT9fi1E7Qt20Hsfe7ObCn3Id6Vg4cr8pjYNDhw6N/Juauh1o0e3eASimpu4EYGrqLgC63bsNfQ8Anc7dhr7HuT81VVxPT99l6J2AYnr6DqDF1NQ7gQ7T0+8Appia2g/MMDV1C0ptodu9kSzbRrt9La3WpbTbV9BuX0WW7aTbvRaYpdu9EaU20e3eAmyg290PdJmaegfQLusyNRWuy9SU5d2lVR2KeMu3FNen7roLDZy+6y50lnH69tvR7TZn3/1udLfL2Xe8g3xqirO33spwZoalt76VfOtWlm+4geH27fT37mW4Zw+Dyy+nf/XVDC+5hOXrryefnWXpppvIN21i8dZbyaenOfv2t5N3u5x95zvR7TZnbr+9KPPOO9FKlbycvNvwZujJ9763oPb6PUU7nLb07rvRwInbb0crVeTXanHmjjuKurzjHeTdLotvfzv5hg0FP5s2sfzmt6DULJ3O9WTZDtrtq2i3L6fVupROZy9Ztp1O580otZWpqbeatnybadueaZd3mfa43bTHHU47VLJV0Fbr9pGytWGDla27jGzdafJ/N9AxcjDF9PTbgQ1GtjbT7d5Els3Rbl9Hq7XL1OUqWq1ddDrXkmXb6HZvMHFvMWn3A1N0u7eZvG8HMs6cmW3sU01+edVKpiu9Xu8LwBfm5+cfNV667rW+QGP3KZRGNE0dn/3sZ/XnPve5Rj4OHnydXbsuReshSrXWlFJ+pSujWkmHKE7YoUMH2blzVy1OLvLKV5VX45RyNLVOMbOMQwcPsnPXLsyi3qXgh9WpycM40hyv7AnpoSNH2Ll9+1h1sXzpMdrFp+feHm+8cYQdO7avkmyNkrE6RtehXhfbT8bB/fffz+c///ng9GWtHQvH7o/jjHhs2P3054tWB4nkCj20Yq+mnMVU0d6Tea0yj/KLZU00y6qlSJGRS0Nhlk5a1lrT0n7S3C4+XZv2mFy2mmSsjtF1kHVZKVa6PGlyLBy6H0yTkJCwPrAipdHkWDh0f0SahISEdYA1dSw84v6qORZOSEg4v/gp2hGakJBwPpCURkJCwkRISiMhIWEiJKWRkJAwEZLSSEhImAhJaSQkJEyEpDQSEhImQlIaCQkJEyEpjYSEhImQlEZCQsJESEojISFhIiSlkZCQMBFWrDR6vd7He73ePb1e7zOR+/eZvwdqYQ/YeystPyEh4fxirR0L3wN825xq3WuuAe7r9XrPU/gITUhIWEdY6dH4TwJ/ZH5bJ8F1/xh7zd+vmfvW+c6n5+fnvzxOAU3+CgEWFhYa41xoXOw8Jv5Wjoudx9Xib6VKY6STYOE3Yz+F1y6oZh2N3z3ZuXPnWIyMG+9C4mLnMfG3clzsPK4Gf+fFEGqWLaXXrvn5+V81S5rttSVLQkLCOkDjTCNirHzB2jEYz0nwPdbjuMnvqFmeHKFasiQkJKwDNCqNBtd8DwI989txLDw/P79gft9X+xzjPcA8lQF0H/CFc2M9ISHhQmBNHQub8Ad6vd7zvV7vWC3NJ8z3XJ9PjoUTEtYX1tSxsFnCzI2TJiEhYX0g7QhNSEiYCElpJCQkTISkNBISEiZCUhoJCQkTISmNhISEiZCURkJCwkRISiMhIWEiJKWRkJAwEZLSSEhImAhJaSQkJEyEpDQSEhImwvnwEer5A21Kk5CQcPFiTX2EGjj+QMdMk5CQcJFipTONT1I44oHKR6jEp+fn5/cZJTFumoSEhIsUa+oj1ED6Ax0nTYn7779/hSwmJCSsJlbsT6MJNa9d907qD/Tzn/+8WhuuEhISzhVr6iM04g90XL+iCQkJFyHW2kdoyB/ofChNQkLC+sCa+ggN+QMdkSYhIWEdQGmtLzQPCecZRokv0PCxql6v95mmj1klXLzo9Xr7Y4PyuDIQwpobQlcbTZVdycM4T/xZG9E++y2Y84n6Ppler7c3JlhmJngvcDE+w/2Y7+WM+3nP1cQEMrj3QjnRNu33BQqzgLw3lgzEsK62kY/xwekLunFsBR/EPp+4qPfJjNmG/9woi70XYRvvp3pR8MKF2rxoy4/cXpEMrCulQXNlL3SHaCp/by2s/kHs84nGfTJm5LlQBuqRz9CM4n8D5ec9z7dNbBwZe8DQvRepzW6ivVIS601pNFV2RQ9jFdD4QezadHU/xZukixHbmqOsGZra8DaKbwDvv0Bnl5ra+FGKGcYxEe+nButNafxUQH4Q+zyjaW/NhZxljIsjtbdwH7/QzNTR6/VmKZ7xvwN+vdfrXYzfKl7RXqn1pjSaKnuhN45N/EHsC4AHqZZFzt4aG2ZOId8HbLsAa/KmZ3iEaq2+QDHzOJ9o4u8+4N8ZA+mngYtGqdXaOCgD42K9KY0mgV/Rw1gFNPEX+iD2ecUYe2u+XHsjMRvIYq3R9Ay/XLs/i7FvnEc0trGFeY4LMvx8wMzAemImVt8/dc57pdbdPg0zAr5A7XVWr9d7xH4/NnT/YuHPNNKXKNa624BfWAdLgfOOMdv4KHDbBXpt3cTfZ8z9bT+N3y1ed0ojISHhwmK9LU8SEhIuMJLSSEhImAhJaSQkJEyEpDQSEhImQlIaCQkJEyEpjYSEhImQlEZCQsJE+P8BVPcOzaNsDKAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set model and likelihood into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initialize figiure an axis\n",
    "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "n = 100\n",
    "test_x = torch.zeros(int(pow(n, 2)), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        test_x[i * n + j][0] = float(i) / (n-1)\n",
    "        test_x[i * n + j][1] = float(j) / (n-1)\n",
    "        \n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    predictions = likelihood(model(test_x))\n",
    "\n",
    "# prob<0.5 --> label -1 // prob>0.5 --> label 1\n",
    "pred_labels = predictions.mean.ge(0.5).float().numpy()\n",
    "# Colors = yellow for 1, red for -1\n",
    "color = []\n",
    "for i in range(len(pred_labels)):\n",
    "    if pred_labels[i] == 1:\n",
    "        color.append('y')\n",
    "    else:\n",
    "        color.append('r')\n",
    "        \n",
    "# Plot data a scatter plot\n",
    "ax.scatter(test_x[:, 0], test_x[:, 1], color=color, s=1)\n",
    "ax.set_ylim([-0.5, 1.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
